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Abstract. We investigate one-dimensional driven diffusive systems where particles 
may also be created and annihilated in the bulk with sufficiently small rate. In an 
open geometry, i.e., coupled to particle reservoirs at the two ends, these systems can 
exhibit ergodicity breaking in the thermodynamic limit. The triggering mechanism 
is the random motion of a shock in an effective potential. Based on this physical 
picture we provide a simple condition for the existence of a non-ergodic phase in the 
phase diagram of such systems. In the thermodynamic limit this phase exhibits two or 
more stationary states. However, for finite systems transitions between these states are 
possible. It is shown that the mean lifetime of such a metastable state is exponentially 
large in system-size. As an example the ASEP with the A0A f± AAA reaction kinetics 
is analyzed in detail. We present a detailed discussion of the phase diagram of this 
particular model which indeed exhibits a phase with broken ergodicity. We measure 
the lifetime of the metastable states with a Monte Carlo simulation in order to confirm 
our analytical findings. 



PACS numbers: 82.40.-g, 02.50.Ga, 05.70.Ln, 64.60.Ht, 64.60.My 
1. Introduction 

The closely related questions of phase coexistence and ergodicity breaking in noisy one- 
dimensional systems are intriguing and have received wide attention in the context of 
driven diffusive systems [DElEllllElEllZlElin!. 

One-dimensional systems with short-range interaction and finite local state space do 
not show phase transition in thermal equilibrium. However, systems far from equilibrium 
may exhibit phase coexistence as well as ergodicity breaking [TJ El El El U\ ■ 

In [10! a new type of phase coexistence was shown to occur in a driven diffusive 
system with weak non-conservation which has no analogue in conserving systems. 
Subsequently it was demonstrated [11] that similar systems may show diverse behavior: 
namely, a slight modification in the transition rates can lead to the appearance of a 
novel phase in which ergodicity is broken. In this paper we give a detailed description 
of a broad class of models including that of PU] and JT] and present a simple criterion 
which can be used to determine whether a given model in this class (described later) 
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shows ergodicity breaking or not. As an example we study analytically and numerically 
the model of ^1] in detail. 

The particular model considered in is the Totally Asymmetric Simple Exclusion 
Process (TASEP) with open boundaries augmented with a creation-annihilation 
dynamics on each lattice site. The TASEP is a well-known and probably the simplest 
nontrivial model of particles with biased diffusion on a one-dimensional lattice. Each 
lattice site can be either occupied (A) or empty (0), double occupancy is not allowed. 
Particles can only hop to the rights with a given rate provided that the target site is 
empty. Particles are injected at the first and removed at the last site at given rates. 

In addition to the biased diffusion of particles on an open lattice the model of ^H! 
includes a simple creation- annihilation dynamics in the bulk 

0^A (1) 

with rates uo a (creation) and Ud (annihilation) on each lattice-site. The mechanism 
is inspired by motor proteins moving along actin filaments, where attachment and 
detachment is possible. Earlier this model was introduced as a toy model reproducing 
stylized facts in limit order markets [T2] . 

The phase diagram of the above system exhibits a phase where a domain of low 
particle density coexists with a domain of high density. The domain wall, i.e. a sharp 
interface between a high and a low density region, is localized in the bulk [HHIT31IT^l ll5|. 
In this phase the two different domains do not represent two possible global steady states. 
The process is ergodic even in the thermodynamic limit and no hysteresis is possible. 
The existence of such a phase was very recently confirmed experimentally in a system 
of single-headed molecular motors (KIF1A) moving along a microtubule jTE] . 

In the following we show that in similar models, i.e. models with biased diffusion 
plus non-conserving reaction-kinetics, an additional phase (or phases) can be built up 
in which ergodicity is broken [TI]. Here the high- and low-density phases do not coexist 
(like in PUj): either the one or the other is present. For finite systems transitions between 
these states are possible but the transition times are exponentially large in system-size. 
The main dynamical mode, which determines the behaviour, is the random walk of a 
shock in an effective potential [llj . This one-particle collective mode provides a unified 
framework for understanding phase coexistence and ergodicity breaking in a certain 
class of systems. Based on this physical picture one can deduce the phase diagram of 
the particular model under consideration. 

Ergodicity breaking is often associated with spontaneous symmetry breaking (SSB). 
However, the former is more general and one can speak about SSB only if the two (or 
more) stationary states are related by some symmetry. In the class of systems we 
consider this is usually not the case (another example for ergodicity breaking but not 
SSB is the Ising model with nonzero external field below the critical temperature). 

The article is organized as follows: In section El we give some introduction to the 
hydrodynamic description we use for a class of reaction-diffusion systems. In section |3] 

| in a more general version (ASEP) particles can also hop in the opposite direction with smaller rate 
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we present the effective dynamics, which is a simple random walk of the shock position 
in an energy-landscape. This physical picture leads to a criterion for the existence of a 
non-ergodic phase. A particular model, which satisfies this criterion is analyzed in detail 
in section |3J Here we also compare the analytical predictions to Monte Carlo results. 
Finally, conclusions are drawn in section ED 

2. Hydrodynamic description 

In the following we restrict ourselves to one-species lattice gas systems in one dimension, 
but we note that considering models with continuous space would not make a big change. 
The description we present here is based on a hydrodynamic approach within Eulerian 
scaling. In this scale the lattice spacing a — > leading to a continuous space-variable 
x = ka and simultaneously the time is also rescaled as t = at m i CTOSCOp i C . Conserving driven 
particle systems on Euler-scale can usually be exactly described by partial differential 
equations for the local particle density jTTj: 



where j(p) is the macroscopic current which one would observe in an infinite system (or 
on a ring) in the stationary state with uniform particle density p. In a continuum 
description of interacting particle systems normally there is a term (referred to as 
diffusive or viscosity term) including a second derivative of the density. Note that 
on Euler scale this term vanishes since it has a prefactor proportional to the lattice 
spacing a. 

It is usually a highly non-trivial question how to determine j(p) for a given model 
and we have to emphasize that the mean-field approximation usually gives a wrong 
answer. In most cases neglecting short range correlations would not reproduce the 
correct j(p) even qualitatively since these correlations play the major role in this 
question. However, there are systems where the homogeneous product measure (i.e. 
a factorized state without correlations) is stationary: in this case the mean-field, by 
coincidence, would give the right answer but this is rather the exception than the rule 
(the ASEP falls into this class). As an analogy we mention that the first order Taylor 
expansion also turns out to be exact when "approximating" linear functions. In this 
analogy improved mean-fields are similar to higher order Taylor expansions. 

Unfortunately there is a wide-spread misunderstanding about equation (J2J) 
according to which this is often referred to as a kind of "continuous mean-field" 
description. The fact that there is a closed equation for the density does not mean that 
correlations are neglected, as one might think, but rather they are built in j(p) which 
is indeed nothing else but a certain correlation function. Another common objection 
is that noise is neglected. Of course, noise is present on microscopic scale but this is 
scaled away on Euler-scale. To summarize: if j(p) is exact then (J2J) is also exact within 
Eulerian scaling (and has nothing to do with any kind of mean-field). 

The systems we consider in this article can be constructed in general as follows: 



dp(x,t) dj(p(x,t)) 
dt dx 



(2) 
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(i) Take a one-species driven diffusive system in one dimension (which can be described 
by (J2J) and couple it to particle reservoirs at the two ends. 

(ii) Introduce non- conserving dynamics in the bulk with reaction rates proportional 
to the inverse system-size. This rescaling of reaction rates is essential to obtain 
nontrivial new behaviour. If these rates are of larger order than 1/L then the 
steady state of the system is governed only by these reactions in the case of L — > oo 
and the boundaries do not play any important role. In the opposite case, when 
they are smaller than (9 (1/L), the system behaves like the original model without 
bulk non-conservation ^HIEIE] (except on the transition line 

As it was argued in ^3] these systems on the Euler-scale (which now include the rescaling 
of bulk reaction rates according to point (|iij)) can be described by the hydrodynamic 
equation 

dp(x,t) dj(p(x,t)) 

~^T + dx = S W*M ® 
where S(p) is some kind of source term coming from the reaction kinetics in the bulk. 
The assumption this description is based on is that in the limit where L —>■ oo the 
probability of a non-conserving process to happen on a finite segment of the chain goes 
to zero. Consequently the system has enough time to relax (locally) to the stationary 
structure of the state which one would obtain with homogeneous density p and without 
non- conserving rates, meaning that the local correlations are the same as those in the 
corresponding conserving system. This implies that the j(p) appearing in is not 
effected by the bulk reaction-kinetics. 

For the correct evaluation of S(p) it is also important to know the local correlations. 
Suppose that in addition to a given driven diffusive particle-conserving dynamics (point 
(0) above) the non- conserving processes from the local configuration C to C with rates 
r(C',C) are included (point (jn))). In the hydrodynamic limit we keep 

R{C',C) =r{C',C)L (4) 

fixed while the system-size L is increased. Let P(C,p) denote the probability of finding 
locally a specific configuration C in the stationary state of the conserving system 
with homogeneous density p. Knowing P(C,p) (which requires the knowledge of local 
correlations) one can easily obtain S(p) as follows: 

S(p) = Y,HC , X)R(C',C)P(C,p), (5) 

c,c 

where k(C',C) simply denotes the particle number difference between configurations C 
and C. 

As an example we give these quantities for the TASEP equipped with the simple 
creation-annihilation dynamics which was considered in |12L ITU1 IT4*] . The TASEP 
falls into the class of models having product measure stationary states (i.e., without 
correlations) on an infinite chain, which implies j(p) = p(l —p). On the other hand, with 
the notation Q a = R(A, 0) and Q d = R(0, A), the source term is S(p) = VL a {l— p) — VL d p 1 
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since P(A,p) = 1 — P(0,p) = p. We note that (a) including partial asymmetry in the 
hopping dynamics would only introduce a prefactor in the current (b) Due to the reasons 
described above in this specific case a mean-field ansatz would also lead to equation (jHJ) 
with the same j(p) and S(p) [TP] , 

Stationary solution 

By requiring time-independent solutions, equation (jH)) reduces to the ordinary first order 
differential equation 

^J'(P(>)) = v c {p{x))p'(x) = S{p(x)), (6) 

where v c (p) = is the collective velocity. This equation allows for the calculation 
of the steady state density profile p(x). As boundary conditions the two reservoirs 
fix the densities at x = and x — 1 (in the hydrodynamic limit we choose a = 1/L 
so the continuous space variable x runs from to 1) but as the differential equation 
(JHJ) is only of first order the solution usually fits at most one of the boundaries. This 
discrepancy is resolved by the appearance of discontinuities which can be located at the 
boundaries (x = 0, 1) or in the bulk. Usually the way of solving (JHJ) numerically with 
given boundary conditions is the following. First one has to regularize the equation 
introducing a viscosity term proportional to e which contains a second order derivative. 
This allows one to solve the resulting equation (since it became second order) and fit 
both boundaries. Finally one lets e — > while usually one or more discontinuities appear 
in the solution. Nevertheless, the resulting curve is a solution of (0) apart from these 
discontinuities. 

However, in practice one does not have to perform the above program for all the 
possible choices of boundary conditions. In [T3] some rules are given to decide which 
are the stable discontinuities at the boundaries. These are the consequences of the 
assumption that the hydrodynamic description is based on, and for conserving systems 
they reduce to the current minimization/maximization principle derived in [TH|. fT^| . 

On the other hand, to have a localized shock in the bulk at xq it is necessary to 
satisfy 

j(p(x -0))=j(p(xo + 0)), (7) 

which follows simply from mass conservation. Additionally, the structure of this shock 
has to be microscopically stable (e.g. a "downward shock" in the ASEP is not stable. 
For details see (THj) 

Given these rules it is still nontrivial to deduce the stationary density profile in 
general and the question remains whether there are cases when there is no unique 
solution, i.e. there are several possible stationary profiles which would imply ergodicity 
breaking. In the following we address this question by studying the dynamics of shocks 
in these systems. 
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3. Shock dynamics — Criterion for ergodicity breaking 

3.1. Site dependent hopping rates 

Suppose that there is a shock in the system which connects two solutions of ©: Picftt^) 
and p r ight(^)- Even if it satisfies (jZJ) and it is a "true" localized shock (note that ((Zj) is 
only a necessary condition) its position fluctuates in time on the lattice scale. According 
to [20] we assume that the shock position performs a simple random walk with effective 
hopping rates Aeft and -D r i g ht which characterize the hopping to the left and right by 
one lattice site and are functions of the densities p r i g ht an d pi e ft - In some special cases 
this effective one-particle description is exact |2U 122 ESj • The average shock velocity 

V s = -Dright — Aeft (8) 

is easy to obtain from mass conservation (the non- conserving processes have different 
time-scale due to the rescaling of rates (JIJ) so they do not enter here): 

v ^ = j(Pright(x)) ~ j(PlcftQr)) , g v 
Pright (%) - Pleft(^) 

and the additional knowledge about the diffusion coefficient D = (-D rigri t + -Dieft)/2 
uniquely defines the effective hopping rates. These rates for the ASEP are 

j(pMt(x)) 



AeftW 



Pright (a;) - Piett(x) 



Dright(x) = i^^L (10) 
Pright W - PleftW 

In a conserving system where the stationary solutions of (|2J) are p(x) =const. the hopping 
rates are also constant all over the system. But as an effect of the nonconservative 
dynamics a source term S(p) appears in (j3J) and leads to non-constant stationary 
solutions, therefore the hopping rates for the shock position also become space- 
dependent. Similar site dependent rates were used in [HI . 
Suppose that there is a point x where v s (x ) = 0, i.e., 

j(pMt(x )) = j(pright(zo)). ( U ) 

There are two possibilities: the shock position at xq can be either stable or unstable. 

(i) Having a stable shock position requires ■^v s (x)\ x=Xo < 0. In this case the shock 
position always has a drift towards the stationary point xq. Differentiating the rhs. 
of (JHJ) and using (0) and (fTT|) the above condition yields 

ff(Pright(so)) ~ S(p M t(x )) < Q ^ 

Pright 

(X ) - Pleft(^o) 

This is equivalent to requiring 

S(p a ) > S( Pb ), (13) 
where p a = min(p right (x ), pieft(^o)) and p b = max(p ri g ht (x ), Pi e ft(^o))- 
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(ii) On the other hand 

S(p a ) < S( Pb ) (14) 

implies an unstable shock position. In this case small deviations from x are 
amplified. 

In the following we show that unstable (fixed) shock positions lead to ergodicity breaking 
in the thermodynamic limit. 

3.2. Effective potential 

The simple random walk (i.e. only nearest neighbour jumps are allowed) of a single 
particle, i.e., the shock position in a confined volume in one dimension necessarily results 
in an equilibrium dynamics which can be characterised by an effective potential. This 
potential on the lattice scale satisfies 

exp(E k - E k+1 ) = — — , (15) 

Qk+l 

where p k is the hopping rate from site k to k + 1, while q k is the one from site k to k — 1. 
It can be shown that the potential in the hydrodynamic limit using the rates Di e { t (x) 
and A-ight i x ) becomes 

E ^^f^{W)) d "- (16) 

In the case of a stable shock position (case (0)) E(x) has a local minimum at Xq. 
The probability density to find the shock at position x is proportional to exp(— E(x)) 
which can be approximated by a Gaussian around Xq. The width of the distribution is 
proportional to L 1 / 2 on the lattice scale (and to L~ x l 2 on Euler-scale) j!4j . 

In the opposite case (case (juj)), when the shock position is unstable, the potential 
E[x) has a local maximum at xq resulting in (at least) two local minima. Since 
E(x) is proportional to L the potential barrier between these metastable positions is 
proportional to L too. From this one expects transition times which are exponentially 
large in L leading to ergodicity breaking in the thermodynamic limit. In section l4~21 it 
is shown that for systems having two metastable shock positions at the two boundaries 
the average transition time between these states grows exponentially with the system 
size (with an algebraic prefactor). It is plausible to assume that the same is true when 
the metastable states are in the bulk. 

3. 3. Criterion for ergodicity breaking 

To conclude, the conditions which has to be satisfied for non-ergodic behavior are 
summarized below: 

If one can find p a and p b such that 

(i) Pa < Pb 

(ii) the shock connecting p a and p b (either p a \pb or Pb\p a ) is stable 
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(iii) according to (JTTJ) j(p a ) = j(Pb) and 

(iv) S(p a ) < 3 (fib), 

then non-ergodic behavior is possible. Otherwise this kind of mechanism does not give 
rise to ergodicity breaking. 

3. 4- Remarks 

To check whether the above criterion is met the specific form of the current density 
relation of the driven diffusive system has to be known. In this subsection we give some 
examples where conditions HHTvl can or cannot be satisfied. 

1. Assume that the conserving part of the dynamics is the ASEP. The current is 
j(p) ~ p(l — p) therefore according to condition lull p& = 1 — p a . The ASEP has a 
product measure stationary state on an infinite lattice (or with periodic boundaries). 
As it is argued above, weak non-conservation with rates of order 1/L (and the open 
geometry) does not change this property in the bulk in the L —>■ oo limit. Because of 
this factorisation every kind of reaction kinetics give a polynomial contribution to S(p) 
the order of which is the number of lattice sites involved in the reaction. It is easy to 
show that with the physical constraint 5(0) > and S(l) < first and second order 
polynomials cannot satisfy condition together withQ One needs to include at least a 
third order source term, and correspondingly a three-site reaction, for having ergodicity 
breaking (in the ASEP condition EH is always satisfied). 

2. It is also easy to see that every driven diffusive system equipped with a simple 
creation-annihilation dynamics (0) is ergodic. The reason is that (pj always gives a 
linearly decreasing S(p) which cannot satisfy condition ITvl 

3. In general, however, it is not necessary to include three-site reactions to get 
ergodicity breaking. As an example it can be shown that in the fc-step exclusion process 
[21] with the two-site reaction 00 ^ A0, the above conditions can be satisfied. In 
this system the steady state is factorized, like in the ASEP, but the current is not a 
symmetric function of the density, therefore the jumps do not have to be symmetric 
with respect to p = 1/2. 

In the following, because of its simplicity, we study in detail the ASEP with the 
reaction A0A ^ AAA. 

4. An example: the model of ref. [llj 

The model we consider in this section was introduced in pi] . It was shown that the phase 
diagram exhibits a phase where ergodicity is broken and the underlying mechanism, 
i.e., the random motion of a shock in an effective potential was also presented. Here 
we present a more detailed discussion of the phase diagram of this model and measure 
the lifetime of the metastable states in the non-ergodic phase to be compared to the 
analytical prediction. 



Ergodicity breaking in one- dimensional reaction- diffusion systems 



9 




Figure 1. The source term S(p) compared to S(l 
There exists a p a < 1/2 such that S(p a ) < S(l — p a ) 



p) for n a = 0.7 and n d = 0.2: 



In the model we consider in the following the driven diffusive dynamics is defined 
by the TASEP, but including backwards hoppings with smaller rate (ASEP) would not 
change the behaviour. The current which enters the hydrodynamic equation is 

j(p)=p(l-p), (17) 

which gives the collective velocity 

v c (p) = l-2p. (18) 

As argued above at least three point reactions have to be coupled to the ASEP to allow 
for unstable shock positions. In what follows we consider the reaction kinetics [TT] 

A0A ^ AAA, (19) 

with rates u a and u^. This may be interpreted as activated Langmuir kinetics, and 
yields a source term 

s(p) = {i a p 2 (i-p)-n dP \ (20) 

where Q aj d = u a ^L are the rescaled rates according to (HJ) . As illustrated in figure [T] for 
certain choice of rates this source term fulfills the conditions for the existence of unstable 
shock positions. At the first lattice site particles are created at rate p_ when this site is 
empty, while at the last site they are annihilated at rate 1 — p+ . These boundary rates 
represent particle reservoirs at the two boundaries with densities p_ and p + . 

In order to calculate the stationary density profiles we rewrite the hydrodynamic 
equation (JHJ) as 

S(P) 



p'(x) 



Vc{p) ' 



(21) 



Integration directly yields the implicit formula for the flow-field (the set of p(x) curves 
satisfying the ordinary first order differential equation ([21))): 



x{p) 



tt a p 



In 



1 

P 



+ c 



(22) 
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Figure 2. The flow-field p(x) of the differential equation J5TJ- The parameters are 
il Q = 0.7 and Qd = 0.1. The solutions are monotonically increasing and convex in 
the regions < p < 1/2 and p* < p < 1 while they are decreasing and concave for 
1/2 < p < p* 1 where p* = Q, a /{Q. a + fid) = 0.875 is the stationary density of the 
reaction-kinetics ()19(l without biased diffusion. Since the collective velocity v c {p) is 
zero at p= 1/2 the solutions are singular at this point. 

where p* = (fl a + Qd)/Q a is the stationary density of the reaction kinetics (fTTJj) without 
biased diffusion and c is an integration constant. The corresponding flow-field p(x) is 
shown in figure |21 The direct calculation of the inverse function is possible in terms of 
the Lambert-functions but unfavorable. 

4-1. Phase diagram 

The solutions (j22j) enable us to deduce the steady-state phase diagram of the system in 
terms of the boundary densities p_ and p + (see figure EJ). 

Before the discussion of the phase diagram let us define some special values of 
the left boundary density, which correspond to special lines of the flow-field. If the 
right boundary density is set to 1/2 there are two solutions which satisfy this boundary 
condition due to the singularity at p = 1/2. We call pi( 2 ) the value of the upper (lower) 
solution at x = (see figure |2J). In addition we define po as the starting point of the 
solution fitting the right boundary density p + = 1 (see figure |2J)- 

First we consider the case p_ < p 2 and 1/2 < p + . In this case pi e ft(^) and p r i g ht(^) 
are the solutions (|2^j) with the boundary conditions pi e ft(0) = p_ and p r i g ht(l) = P+ 
and they are well-defined in the whole range of x. As Pnght^) > 1/2 > pi e ft(^) the 
shock connecting the two solutions is stable so one can define the effective potential 
E(x) for the shock position according to (fTp^l . In figure H] four qualitatively different 
energy-landscapes are shown which are present in the region of the phase diagram 
under consideration. They are calculated by inserting the analytical solutions of pi e ft(^) 
and pieft(a;) in into (fl"6|) . These effective potentials indicate four qualitatively different 
behaviours: 
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Figure 3. Phase diagram for fixed fl a — 0.7 and S! f j = 0.1 with two high-density 
phases (HDI, HD2), a low-density phase (LD), a coexistence phase, and the non- 
ergodic phase. 
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Figure 4. Examples of the effective potential in four phases. Note that in the HD 
and LD E(x) can be either convex or concave. Parameters are: il Q = 0.7 fid — 0.1; 
P- = 0.25 p + = 0.6 (non-ergodic) ,p- = 0.1 p + = 0.903 (coexistence), p- = 0.27 p + = 
0.6 (HD), p_ = 0.1 p+ = 0.89 (LD). 




• There is a region in the parameter space where E(x) is monotonically decreasing 
resulting in a stationary shock position at the right boundary. Consequently the 
bulk solution is given by pi e ft(x) - this is the low density phase (LD): 

• In another region, called high density phase (HD), the potential E(x) is 
monotonically increasing, therefore the stationary shock position is at the left 
boundary and the solution is Pright^) in the bulk. 

• In the coexistence phase E(x) has a local minimum in the bulk leading to a 
stable shock position at a macroscopic distance from both boundaries. The shock 
separates a low density - and a high density region characterised by pi e ft(^) and 
Pright(^)- This phase is analogous to the one found in [10]. When measuring the 



Ergodicity breaking in one- dimensional reaction- diffusion systems 



12 



density profile by a time average one observes a smooth curve connecting the two 
regions. This is the consequence of the fluctuation of the shock position, the shock 
itself is sharp. The width of this "transition region" can be calculated by a harmonic 
approximation of the potential in the vicinity of the minimum, which gives an 

£1/2 

dependence on the lattice scale (on the Euler-scale the width vanishes as L -1 / 2 ). 

• In addition a "non-ergodic" phase is found, which is characterised by a potential 
with a local maximum in the bulk. The two minima at the left and right boundary 
correspond to the two meta-stable states: pi e ft(^) and p T i g ht(x). For finite systems 
a transition between the low-density and the high-density state is possible but very 
improbable if L is large due to the potential barrier AE ~ L (see section l3~2*)) . Of 
course strictly speaking for any finite L there is an unique stationary state which 
is some kind of linear combination of the high - and low density states but the 
relaxation time to this stationary state diverges exponentially with L. In the limit 
L —>■ oo the system becomes non-ergodic since one of the stable states is selected in 
a very short while (depending mainly on the initial configuration) and a transition 
to the other one is not possible, therefore a time-average does not reproduce the 
ensemble- aver age. 

With the picture of the moving shock in mind it is also possible to derive exactly 
the two lines in the phase diagram shown in figure El which separate these phases. One 
of them is determined by the condition E'(0) = 0. To the right (left) of this line a shock 
at the left boundary - and therefore the high density state - is stable (unstable). Along 
the other line E'(l) = 0, and it separates the two regions where the shock at the right 
boundary - and therefore the low density state - is stable (to the left) or unstable (to 
the right). The coexistence phase is the one where neither of them are stable and as a 
consequence the shock is driven into the bulk and is localised there. In the non-ergodic 
phase both the low density and the high density state are stable. 

The sign of the slope of the effective potential is determined by the sign of the 
average shock velocity v s (x) (see equation (jSESJ) so in what follows we focus on this 
quantity at the two boundaries. The shock velocity at the left boundary v s (0) is zero if 
j(pieft(0)) = j(Pright(0)) and therefore 



where PrightC^) is the solution fitting the right boundary density p + . Increasing the 
value of p_, v s (0) decreases so the high density phase is stable to the right of the line 
defined by (J2HJ)- This phase boundary starts at p_ = 1 — p ,p+ = 1 and ends at 
P- = l-pi,p+ = 1/2 

The analogous equation for the phase boundary of the stable low-density state is 



where PMt(x) is the solution fitting the left boundary density p_. For similar reasons 
the low-density state is stable to the left of this line, which starts at p_ = 0, p + = 1 and 



P- = Plcft(O) = 1 - Pright(O) 



(23) 



Pleft(l) = 1 - Pright(l) = 1 - P+, 



(24) 
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Figure 5. Shown is the narrow phase near the intersection point of the two nontrivial 
curves of the phase diagram (figure 0J- The two thick solid lines represent these phase 
boundaries in this schematic enlarged figure. Along the thin solid lines the energy 
profile has a stationary point at a given x (x\ > X2 > £3). In the coexistence phase 
this is a local minimum whereas in the non-ergodic phase this is a local maximum. 
The second derivative changes sign on the thick dashed curve. In the "triangular" area 
the energy landscape has a local minimum and a local maximum as well in the bulk. 

ends at p_ = P2,p+ — 1/2. Note that the condition for having a non-ergodic phase is 
p 2 > 1 — Pi, which is not satisfied for all values of Q a ,d- 

We note that the analytical treatment suggests an additional phase in the vicinity of 
the intersection point of the two lines discussed above (see figure EJ). In this phase E(x) 
has a local maximum and a local minimum as well (see figure EJ). This type of energy- 
landscape also leads to ergodicity breaking. In one of the stationary states the shock is 
located in the bulk, in the other it is at the right boundary (LD). Although this phase 
is very narrow and not observable by Monte-Carlo simulations in our specific model, 
it indicates that in general E(x) can be more complex and can have more stationary 
points. 

Now let us consider another region of the phase diagram. Suppose that keeping 
p + fixed between 1/2 and 1, we increase p_ beyond p 2 . Due to the singularity the 
left solution p\ c {t{x), and therefore E(x), is not well-defined in the whole range of x 
anymore (see figure Ej). Nevertheless it is perfectly defined in a macroscopic interval 
x G (0, z), where < z < 1 is the point where p\ e ft{x) becomes singular. However, due 
to the negative average speed the shock stays in this region so there is no change in 
the qualitative behaviour. If p_ approaches 1/2 the length of the interval where E(x) 
is defined shrinks to zero. Further increasing of p_ leads to the loss of the effective 
potential. Although the bulk behaviour does not change in this case we expect a change 
in the details in the way the density profile approaches the bulk solution on the lattice 
scale near the left boundary. 

The situation is similar to the case in the (conserving) ASEP in the high density 
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Figure 6. The potential in the narrow phase near the intersection point of the 
lines separating the HD, LD, coexistence, and non-ergodic phases. Due to the small 
energy-barrier one would need huge systems to observe this phase by a Monte-Carlo 
simulation (compare the scale on the "y" axis to the one of figure @J|. Parameters are: 
fl a = 0.7 fl d = 0.1 p- = 0.138 p+ = 0.8435 



phase. If the left boundary density p_ is less then 1/2 there is a "pure" shock localized 
at the left boundary leading to a simple exponential decay to the bulk density, which is 
equal to the right boundary density p+. However, if the left boundary density is above 
1/2 there is an algebraic prefactor to the exponential l2*o]. which can be interpreted 
as a result of a combination of an algebraic decay to the maximal current state p = 1/2 
and a shock from this value to the bulk p = p + (for details see [20J). Similarly, in our 
model we expect a pure exponential in the approach to the bulk p(x) in the HD for 
p_ < 1/2, while for p_ > 1/2 an algebraic modification is expected. 

Nevertheless the bulk behaviour, which is given by p r i g ht(x), does not change in the 
whole region p + > 1/2, p_ > p 2 so this belongs to the HD. 

The above situation is analogous to the one where being in the LD we decrease p + 
with fixed p_. Although E{x) could be defined also in the case of p + < 1/2 since both 
solutions pieft(^) and p r i g ht(^) are regular, an algebraic modification of an exponential 
decay to the bulk profile is expected, just like in the usual ASEP in the LD if p + < 1/2 
[2*H| EE] ■ The bulk behaviour however, remains unchanged and is governed by p_ in this 
region so this belongs to the LD. 

Finally suppose that the system is in the HD discussed before and we decrease p + . 
As p + approaches 1/2 the bulk profile converges to the solution which starts at pi and 
ends at 1/2. We argue that further decreasing of p + does not change the bulk behaviour. 
This case is analogous to what happens in the conserving ASEP in the HD-maximal 
current (MC) phase boundary. There, an algebraic decay to the bulk density 1/2 sets in 
if p + < 1/2. Here, we expect similar behaviour near the left boundary, but it is still not 
trivial to predict the type of decay because this algebraic decay is somehow combined 
with the singularity of the bulk profile at x — 1. Consequently, in the whole range of 
parameters p + < 1/2 and p_ > p 2 the bulk stationary profile is independent of both 
boundaries (and is given by the thick line starting at pi on figure |2J). In this sense this 
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is similar to the MC phase in the ASEP phase diagram and therefore we distinguish 
between the HD1 (p + > 1/2) and HD2 (p + < 1/2) phases. 

The region 1 — p\ < p_ < p 2 and p + < 1/2 belongs to the non-ergodic phase. The 
stationary state here is either the corresponding low-density state determined by p_ or 
the state corresponding to the HD2. Summarizing, we can say that the bulk behaviour 
below the p + — 1/2 line is just the same as on this line itself. 



4-2. Transition times 

The existence of two metastable states is not a sufficient condition for ergodicity 
breaking. Here we show that the transitions from one state to the other are impossible 
in the thermodynamic limit. This is done by showing that the mean first passage times 
(MFPT) of the shock from one boundary to the other grows exponentially with the 
system size (which is suggested by the linearly growing potential barrier). 

The MFPT r, of a random walker with site-dependent hopping rates from site 1 to 
L can be calculated using a formula derived by Murthy and Kehr |2*7j . 

L-1 L-2 L-1 i 

^ = E 1 + E 1 EII- ( 25 ) 

k=i P k k=i P k i=k+i j=k+i Pi 
where pj is the hopping rate of the random walker from site j to j + 1 and qj from site 
j to j — 1. According to (|T3|) this can be written as 

. „Ei—E k 

r= £ V (26) 

l<fc<Z<L-l ^ K 

This formula reduces only in some special cases to a closed form and one has to rely on 
numerical evaluations in general. We tested the scaling of r with L for several shapes 
of the energy landscape with a maximum proportional to L. We found that the MFPT 
r from site 1 to site L is exponential in L with an algebraic prefactor: 

r ~ L a exp(AE), (27) 

where AE is the maximum of E(x) — E(y) with the condition x > y, and is proportional 
to L. The value of a depends on the specific form of E(x), more precisely on the 
behaviour near the minimum and maximum of it. In our case, when at the minimum 
(at the boundary) E(x) starts linearly and at the maximum it is quadratic a = 1/2. 

The formula (|2*7jl is also supported by an analytical argument based on 
approximating the double sum of by a double integral: 

r f L{E{x)-E(y)) 

r^L 2 / — -— dxdy. (28) 

J0<y<l J0<x<y -bright [U ) 

Here we introduced E(x) = E(x)/L, which remains constant with increasing L and we 
used the right hopping rate A-ight of the shock. These are functions of the continuous 
space- variable x. For energy profiles similar to the solid line of figure 0] the main 
contribution to (|28|) comes from the part where y is close to zero and x is in the vicinity 
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of the local maximum of E(x). Using the first and second order Taylor-expansions of 
E(x) near zero and Xq (the point where E(x) has its local maximum) and extending the 
limits of the integrals one gets 



This is indeed of the form of (J27j) with a = 1/2. We note that for energy profiles which 
have different behaviour near the minimum and maximum, the same argument would 
typically lead to the same form with different value of a. 

The MFPTs can be calculated also numerically, using the hopping rates (jl(J|) which 
result from the density profiles suggested by the hydrodynamic description. The energy 
profiles and the resulting MFPTs for different system sizes are shown in figure First 
the hydrodynamic equation (|21|) was solved numerically with the two different boundary 
conditions: p(0) = p_ and p(l) = p + which gave the two solutions pi e ft(x) and p T i g ht(x). 
After discretisation of x the left and right hopping rates were defined by ()1()J) which were 
used to calculate the energy profile (fT3j) and r (|23j) . It can be seen that the height of 
the energy barrier increases linearly with system size and this results in an exponential 
increase of the first passage times. It is also possible to check the (less important) 
algebraic prefactor. In figure ^rexp(-AE) is shown against L and the prefactor L 1 ^ 2 
is confirmed. However, the constant factor in (fHU|) does not agree with the numerical 
results. It can often happen when a sum is replaced by an integral that the functional 
dependence of a parameter is reproduced correctly but not the constant factor. FiguresEI 
and |H1 confirm that the integral approximation (|2*Hj) , resulting in ()30|) , correctly describes 
the behavior of (|25jl for L —>■ oo up to a constant factor. A comparison of this predicted 
mean first passage time and simulation results of the full system is made in section 14.41 

We note that of course the potential E(x) characterizes only the equilibrium 
behaviour of the random walker: it gives the probability distribution of the shock 
position along x after relaxation, so for times f > r. In contrast, the MFPT is a 
quantity which refers to non-equilibrium and is not given by E(x) only. Note that r 
is a function of the hopping rates which are not uniquely defined by E(x). However, 
it turns out that the behaviour of r for large L pTj) does not depend strongly on the 
specific choice of the hopping rates, but it changes only the prefactor (gives a different 
time-scale) . 

To summarize, we conclude that the transition times of the shock from one 
boundary to the other diverge exponentially with system size and therefore ergodicity 
is broken in the thermodynamic limit. 




(29) 



For large L the y-dependence of -D right can be neglected leading to 



r 




(30) 
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Figure 7. Energy profiles for p_ = 0.2705, p+ = 0.63, n a = 0.5, fl d = 0.1 and 
different system sizes. The inset shows the mean hrst passage times from left to right 
(solid line) and from right to left (dashed line) against the system size. 
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Figure 8. r/ exp(A_E) is shown as a function of L for the LD (x) and HD (+) states. 
The parameters are: tt a = 0.5, fid = 0.1, p- = 0.2705, p+ = 0.63. The straight lines 
have slope 1/2 and are fitted by hand (note the logarithmic scale). This confirms the 
L 1 / 2 prefactor in r, which corresponds to a — 1/2. 

4-3. Simulation technique 

In order to confirm our predictions we performed Monte-Carlo (MC) simulations. In 
the effective dynamics described above the only degree of freedom is the position of the 
shock, which has to be defined on the lattice scale to be able to compare the analytical 
results with simulations. 

In the ASEP this is usually done by introducing second class particles (B) |28j . 
These particles behave as normal (A) particles with respect to holes but they count as 
holes in an environment of A particles. If there is only a single B particle in the system 
the position of it is a good definition of the shock position on the lattice scale. In non- 
conserving systems, however, one has to take care of not losing the second-class particle 
(and not introducing an additional one). To this end we introduce another method to 
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Figure 9. Hopping processes in the bulk. All of the above processes have rate one. On 
the upper lattice the dynamics of the original model is recovered while a single particle 
moving on a separate (lower) lattice was introduced to indicate the position of the 
shock. The reactions (• o • ^ • • •) take place only on the upper lattice independently 
of the configuration of the lower one. Note that in our simulations the last hopping 
process does not happen since we introduce only one particle on the lower lattice. 

trace the shock position on the microscopic scale. 

Consider an additional one dimensional lattice with N sites parallel to the original 
(upper) one. The hopping processes are those explained in figure El while the bulk non- 
conserving part of the dynamics (• o • ^ • • •) as well as the injection and ejection 
at the boundaries act only on the upper lattice. One can see that the dynamics of the 
original model is recovered on the upper lattice. 

Notice that without non-conserving processes in the bulk and identifying 

• o • o , 

-> A, -»• 0, B, -> B 31 

o o • • 

one recovers the ASEP with second-class particles §. The above postulated dynamics 
also serve as an alternative definition of the second class particle in an open system 
(there is a minor difference between these dynamical rules and those of [22J at the 
boundaries). 

Introducing only one single particle on the lower lattice makes it possible to trace 
the shock position without having any effect on the (non-conserving) dynamics of the 
original system. Furthermore the setup described above has another practical advantage, 
namely, it is suitable to adapt multispin coding 29 , which speeds up the simulations 
remarkably. 

4-4- Monte Carlo Results 

In the simulations we measured the time dependence of the position of the shock (the 
particle on the lower lattice) in the broken ergodicity phase. A time window of the 
result is shown in figure E3 One can see that the shock spends most of the time in the 

§ Alternatively one can consider the configuration (°) as a distinct (C) type of particle. Since this 
counts as an empty site if it meets A or B but behaves as a particle in an empty environment C could 
be called a third-class particle. 
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Figure 10. Snapshot of the time evolution of the scaled position of the second class 
particle for L = 1000, p_ = 0.2705, p+ = 0.63, n a = 0.5, O d = 0.1. A position of the 
second class particle near the left boundary (x w 0) corresponds to the high density 
state, while a position near the right boundary (x w 1) corresponds to the low density 
state. 
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Figure 11. Numerically determined distribution of the transition times from the high 
density to the low density state (histogram) compared to the exponential distribution 
(dashed line). Similar results are found for the transition times in the other direction 
but with a different mean. Parameters are: L = 1000, p_ = 0.2705, p + = 0.63, il a = 
0.5, ild = 0.1. The histogram is based on a list of 2094 measured transition times. 

vicinity of either boundary: the shock being at the left/right boundary means that the 
system is in the HD/LD state. 

We also measured the first passage times of the shock from left to right and vice 
versa which correspond to the lifetime of the HD and LD states. The one-particle picture 
implies that the first passage times have to follow an exponential distribution which is 
confirmed by our results (see figure [TTJ). 

By measuring the mean first passage times it is possible to check the predictions of 
section 14.21 However, it can be carried out only for relatively small systems because of 
the exponentially increasing time scale. We measured the mean first passage time up to 
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Figure 12. The measured average transition times compared with the theory. The 
solid/dashed line shows the lifetime of the high-density/low-density state according to 
the theory as a function of the system size, and the squares/triangles are the results 
of the MC simulations. Parameters are: p_ = 0.2705, p+ = 0.63, Sl a = 0.5, = 0.1. 
For each value of L the number of observed transitions was between 1600 and 2100. 

L = 1000 (see figure [EJ. The agreement for the passage time from left to right (this is 
the life time of the HD state) is good, but it seems that the theory gives a wrong result 
for the mean lifetime of the LD state (note, that there is no particle-hole symmetry in 
the system). In the following we present a possible explanation for this discrepancy, 
which we believe to disappear as L — > oo. 

For finite systems the hydrodynamic description is not exact, there are finite 
size effects. Our assumption was that the reaction kinetics do not change the local 
correlations, but only the value of p, which becomes x-dependent on Euler scale. This 
is not true in general for finite systems where the reactions modify the correlations 
and therefore S(p) and j(p) is changed in ©. Moreover, the effective hopping rates 
of the shock can also change if there are correlations. In our model, however, there is 
a special point where the stationary state remains a product measure (the stationary 
distribution of the ASEP) even for finite systems. The homogenous product state with 
density p* = Q a /(Q a + Q d ) is a stationary distribution of the process (fT§j) alone and 
also together with the ASEP hopping dynamics. Hence it is plausible that the further 
we are from p* the more pronounced the finite size effects are (note that p(x) = p* is 
always a solution of (|2Tj) ). 

To confirm this we measured the equal time density- density correlations 

Cfc = (n L / 2 n L / 2+ k) - (n L / 2 }{n L /2+k) (32) 

in the high and low density states (ni is the occupation number of the Ith site). The 
results show that in the HD state c& is very close to zero (it is zero within the statistical 
error of our measurement), whereas in the LD state there is a considerable correlation 
in finite systems. Although the correlation length (on the lattice scale) increases the 
amplitude decreases with increasing system size. 
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Figure 13. (Color online) Shown is the measured density profile of the low density 
state for different system sizes compared to the hydrodynamic limit. The average was 
taken over those configurations of the history where the second class particle was on 
the last (Lth) lattice site. It can be seen that the finite size effects are very pronounced 
even for large systems. For comparison the inset shows the measured density profile 
of the high density state for L = 250 and L = 1000 compared to the hydrodynamic 
limit. Here, apart from the exponential behaviour at the boundary (which is due to 
the microscopic structure of the shock) the agreement is perfect even for relatively 
small systems. Parameters are: p_ = 0.2705, p+ = 0.63, fl a = 0.5, fid = 0.1, which 
gives p* = 5/6 w 0.83. 



This suggests that the density profile of the high density state is much less sensitive 
to finite size effects than the one of the low density state. The deviation from the 
"true" p(x), by assuming no correlations, increases with increasing/decreasing x in the 
LD/HD state since they are determined by the left/right boundary density, therefore 
the agreement with the analytical prediction is the worst in the LD state near the right 
boundary (see figure IT3|) . The dominant contribution to the mean first passage time 
from left to right is determined by that part of the energy landscape which is to the left 
of the maximum. In this region the density profiles, and therefore also E(x), are less 
sensitive to finite size effects hence there is a good agreement. On the other hand the 
energy landscape determined by the hydrodynamic description is much less reliable near 
the right boundary, which results in a significant difference between the predictions and 
the measurement of the life time of the LD state for accessible system-sizes (see figure 
[TjZJ). We believe that these are finite size effects and due to the weakening of correlations 
there is agreement for L — > oo (compare the range of L in figure El and figure ITHJ) . 

We also found hysteresis ^T] by measuring the space averaged density p while 
changing p_ in such a way that the system starting from the LD phase passed through 
the non-ergodic phase and ended up in the HD phase. Then the process of changing p_ 
was reversed. Results are shown in figure El We found that the hysteresis loop inflates 
with increasing speed of sweeping, which is reminiscent of hysteresis in usual magnetic 
systems. 
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Figure 14. Hysteresis plot for L = 2000, Q a = 0.7, £l d = 0.1, p + = 0.45. p- was 
changed by 10~ 4 in every 5000 (solid line), 1500 (dashed line), and 500 (dotted line) 
MC steps (time average was not taken). The hysteresis loop gets wider as the speed 
of changing p_ is increased. 

5. Conclusions 

We found a class of noisy one-dimensional systems with only one species of particles 
which can show ergodicity breaking in the thermodynamic limit. These systems can be 
constructed by taking a conserving one-species driven diffusive system on an open one- 
dimensional lattice and introducing bulk non-conservative reaction-kinetics with rates 
proportional to the inverse system size. A general property of these systems is that the 
stationary density profile is not constant in the bulk, moreover it is often non-continuous, 
i.e., a shock is localized at a macroscopic distance of the boundaries. 

By using a hydrodynamic description on the Euler-scale we find the dominant 
dynamical mode, which is the random walk of a shock in an effective potential (within 
a confined volume). It is surprising that in these non-equilibrium many-body systems 
the major role is played by a single-particle collective mode moving under equilibrium 
conditions. In the above framework localized shocks correspond to local minima of 
the potential, while a local maximum generates ergodicity breaking, since the effective 
energy barrier is proportional to the system size L. Based on this physical picture we 
give a criterion for the existence of a non-ergodic phase. In this phase there are more 
stationary states all of which correspond to a local minimum of the potential. In finite 
systems transitions between these states are possible but the average transition time is 
exponentially large in the system-size. 

We study one particular model in detail which is the ASEP with the A0A ^ AAA 
reaction kinetics jllj . With the moving shock in mind it is possible to derive the phase 
diagram of the system analytically, which indeed exhibits a phase with broken ergodicity. 
We also performed Monte Carlo simulations to confirm our findings and the results are 
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in good qualitative and quantitative agreement with the theory apart from the average 
transition time from the low density to the high density state in the non-ergodic regime. 
We argue that this is due to strong finite size effects. Unfortunately the accessible 
system sizes are very much limited in these simulations because of the exponentially 
growing time scale. An analytical calculation of correlations present in finite systems 
would be needed for a better understanding of this discrepancy. 

Ergodicity breaking is often associated with spontaneous symmetry breaking (SSB) 
but we would like to emphasize that the former is much more general. In the systems 
we consider, usually there is no symmetry breaking since there is no symmetry which 
would relate the different stationary states to each-other. However, in special cases the 
ergodicity breaking we observe does reduce to SSB. One example is the ASEP with the 
reaction kinetics A0A — > AAA and 0A0 — > 000 . If these reactions have the same 
rate Q and the two boundary densities are related by p + = 1 — p_ then the exchange of 
particles and holes followed by a space reflection is a symmetry. According to © the 
source term of this model is 

5(p) = -np(l-p)(l-2p), (33) 

which satisfies the criterion for ergodicity breaking (see sec. I3.2|) . The two solutions of 
(JHJ) fitting the left and right boundary densities are 

Pieft(x) = (l + (y - lj e n ^ , (34) 

p rig ht(x) = (l + (J- - l) t^-A . (35) 

These density profiles are related by the above symmetry and generate an effective 
potential (jl(JH6j) which is symmetric and has a single maximum at x = 1/2. For 
p_ = 1— p+ < 1/2 the theory predicts two stationary states which correspond to the 
shock being at the left and right boundaries and are characterised by the above density 
profiles 



This is a novel kind of mechanism for SSB as compared to the bridge model pQ, 
where the SSB is triggered by blockages at the boundaries. These blockages can be 
considered as effective impurities [SU] and induce the breakdown of the hydrodynamic 
description by inhibiting the identification of boundary densities in a consistent way. 
Once these blockages are removed the system becomes ergodic and no symmetry 
breaking is possible. In contrast, in the systems we consider there are no such 
impurities. The boundary densities are well-defined and the hydrodynamic description 
works. However, in some cases the hydrodynamic equation has several stable stationary 
solutions. The dynamical picture behind is the random motion of a shock in an effective 
potential, which has several local minima. 



|| For p_ = 1 — p + > 1/2 the boundary densities have to be replaced by 1/2 in equations 
because of a boundary layer which is similar to that of the ASEP in the maximal current phase. 
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